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^ ■ Abstract 

We present a non-perturbative computation of the running of the coupling a s in 
QCD with two flavours of dynamical fermions in the Schrodinger functional scheme. 
We improve our previous results by a reliable continuum extrapolation. The A- 
parameter characterizing the high-energy running is related to the value of the cou- 
pling at low energy in the continuum limit. An estimate of Ar is given using 
large-volume data with lattice spacings a from 0.07 fm to 0.1 fm. It translates into 
A^ = 245(16) (16) MeV [ assuming r = 0.5 fm]. The last step still has to be im- 
proved to reduce the uncertainty. 
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1 Introduction 



The QCD sector of the Standard Model of elementary particles constitutes a renor- 
malizable quantum field theory. After determining one mass parameter for each quark 
species and the strong coupling constant (at some reference energy) which fixes the 
interaction strength, QCD is believed - and so-far found - to predict all phenomena 
where only strong interactions are relevant. Due to asymptotic freedom, this pro- 
gramme could be successfully implemented for processes with energies E ^> 1 GeV 
by using perturbation theory to evaluate QCD. A recent report on such determina- 
tions of the coupling implied by a multitude of experiments is found in pQ. 

At smaller energy the lattice formulation together with the tool of numerical 
simulation is used to systematically extract predictions of the theory. Here the free 
parameters are typically determined by a sufficient number of particle masses or 
matrix elements associated with low energies, and then for instance the remaining 
hadron spectrum becomes a prediction. See [2] for a recent reference on such large 
scale simulations. 

While this almost looks like two theories having disjoint domains of applicability, 
we really view QCD as one theory for all scales. Then the adjustable parameters 
mentioned above are not independent but half of them should in fact be redundant 
or, in other words, predictable. The ALPHA Collaboration is pursuing the long-term 
programme of computing the perturbative parameters in the theory matched with 
Nature at hadronic energies. Here a non-perturbative multi-scale problem has to be 
mastered: hadronic and perturbative energies have to be covered and kept remote 
from inevitable infrared and ultraviolet cutoffs. To nevertheless gain good control 
over systematic errors in such a calculation, a special strategy had to be developed; 
it will be reviewed in the next section. 

As done for many other applications of lattice QCD, also our collaboration has 
recently set out to overcome the quenched approximation and include vacuum fluc- 
tuations due to the two most important light quark flavours. In this article we pub- 
lish detailed results on the energy dependence of a non-perturbative coupling from 
hadronic to perturbative energies and thereby connect the two regimes of QCD. This 
connection has been broken into two parts. It is first established with reference to a 
scale L max in the hadronic sector that is not yet directly physical but rather techni- 
cally convenient within our non-perturbative renormalization scheme. This part is, 
however, a universal continuum result and represents what is mainly reported here. 
In a second step, which does not anymore involve large physical scale ratios, one 
number relating our intermediate result to physics, like for instance L max F n , has to 
be computed. Here we can cite in this publication only a first estimate and not yet 
a systematic continuum extrapolation, which is left as a future task. 
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Earlier milestones of our programme consisted of the formulation of the finite- 
size scaling technique [HJ, its adaptation to QCD |4|5| . check of universality jHj and 
the complete numerical execution of all steps in the quenched approximation |7|8|9j . 
A shorter account of the present study with a subset of the data was published in [TU] . 

We would like to mention that beside our finite- volume technique there are also 
efforts to compute the coupling and quark masses in a more direct "one-lattice" 
approach. Some recent references with unquenched results include |lipi2piH| . Since 
it is not possible to accommodate all scales involved on one lattice in a satisfactory 
fashion, in these works a perturbative connection is established between the bare 
coupling associated with the lattice cutoff scale a and a renormalized coupling at 
high energy in a continuum scheme such as MS. In our view, this step represents 
the main weakness of the method. 1 Series in the non-universal bare coupling are 
usually worse behaved than renormalized perturbation theory. There are techniques 
to improve this (tadpole improvement, boosted perturbation theory) which, in cases 
where a non-perturbative check is possible, sometimes help more and sometimes help 
less. Thus it is not easy to verify that one uses this perturbation theory in a regime 
where it is accurate and to quote realistic errors on this step. On the other hand, this 
approach is much easier to carry through and typically yields approximate results 
before an application of our strategy is available. 



A non-perturbative renormalization of QCD addresses the question how the high- 
energy regime, where perturbation theory has been successfully applied in many 
cases, is related to the observed hadrons and their interactions at small energies. This 
relation involves large scale separations and thus is difficult to study by numerical 
simulation. Naively, it would require simulations with a cutoff a~ l much larger than 
the largest energy scale, combined with a large system size L, much larger than the 
Compton wavelength of the pion. In summary this would mean 



to avoid discretization and finite-size errors. This clearly corresponds to - with the 
current computing power - inaccessibly large lattices. 

Our method to overcome these problems has been developed in |3|4|15|16|17| 
I5|1U| : pedagogical introductions can be found in |18|19| . The key concept is an 
intermediate finite-volume renormalization scheme, in which the scale evolution of 

As for ref. an additional relevant question concerns the correctness of the formulation of 
fermions on the lattice. We refer to |14) and references therein for an introduction to the problem. 



2 Strategy 
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the coupling (and the quark masses) can be computed recursively from low to very 
high energies. At sufficiently high energies, the scale evolution is verified to match 
with perturbation theory and there the A-parameter is determined. 



2.1 Renormalization 

Any physical quantity P should be independent of the renormalization scale /i. This 
is expressed by the Callan-Symanzik equation |2U|21j 

4 + ^l +Tft) ^} p = - (2 ' 2) 

Here, the /3-function is given by 

m=fi$- (2.3) 

Hence, once the coupling has been defined non-perturbatively for all scales (see sec- 
tion also the /^-function is defined beyond perturbation theory. 2 For weak cou- 
plings or high energies only, the /^-function can be asymptotically expanded as 

m = -f(bo + b 1 g 2 + b 2 f + ...) . (2.4) 



In the following we only consider mass independent renormalization schemes 
in which the renormalization conditions are imposed at zero quark mass. Particular 
examples are the Schrodinger functional scheme described below, and the MS scheme. 
If in two such schemes the couplings can 3 be mutually expanded as Taylor series of 
each other (once they are small enough), 

g' 2 = g 2 (i + c^g 2 + ...) , (2.5) 

then the 1- and 2-loop coefficients bo and b\ in ()2.4|) are universal, 

b o = 77^ (ll-|iVf J , (2.6) 



(An) 2 V 3 

h = -L-fm-^-Nt) , (2.7) 



(4tt) 4 V 3 



2 For the r-function similar statements and expressions are valid, once the running quark mass 
to has been defined non-perturbatively. 

3 See |23) for an example that this restriction can be non-trivial for non-perturbative couplings. 
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while the higher-order coefficients depend on the choice of the coupling. Starting 
from the 3-loop coefficient in the MS scheme [21] and its conversion to the minimal 
subtraction scheme of lattice regularization |25I26|27|28| , the 3-loop coefficient in the 
Schrodinger functional scheme (with all the parameters as specified below) 

b 2 = (0.483(7) - 0.275(5)iV f + 0.0361(5)iV f 2 - 0.00175(l)iV f 3 ) /(4vr) 3 (2.8) 

could be obtained in |2"§j . 

Now, a special exact solution of the Callan-Symanzik equation ()2.2|) is the renor- 
malization group invariant A-parameter, 



A = ii (bof^y^^e- 1 ^ 2 ^ exp \ - dx 



o 



1 1 bx 

+ 



j3(x) b x 3 b^x 



• (2-9) 



A is scale independent but renormalization scheme dependent. The transition to 
other mass independent schemes is accomplished exactly by a 1-loop calculation. If 
in another scheme the coupling is given through (J2.5|) with both g' and g taken at 
the same //, the A-parameters are related by 

A' = Ae c "'/«. (2.10) 

In particular, the transition from the Schrodinger functional scheme with two dynam- 
ical flavours to the MS scheme of dimensional regularization is provided through t 3Qj, 

A^| = 2.382035(3)A (2) . (2.11) 



2.2 Schrodinger functional 

Since we want to connect the perturbative regime of QCD with the non-perturbative 
hadronic regime, we have to employ a non-perturbative definition of the coupling. 
Furthermore, the definition of the coupling should be practical. This means that one 
has to be able to evaluate it on the lattice with a small error, that cutoff effects are 
reasonably small and that its perturbative expansion to 2-loop order is computable 
with a reasonable effort. In principle there is a large freedom for the choice of the 
coupling, however, it turns out that the conditions above are hard to fulfill simulta- 
neously. 

To this end we consider the Schrodinger functional, which is the propagation 
amplitude for going from some field configuration at the time xq = to another field 
configuration at the time x$ = T. Here the space-time is a hyper-cubic Euclidean 
lattice with discretization length a and volume Txll We choose T = L so that 
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L is the only remaining external scale in the continuum limit a — > of the massless 
theory. 

The SU(3) gauge fields U are defined on the links of the lattice, while the fermion 
fields are defined on the lattice sites. The partition function of this system is given 
by 

Z = e- r =[ D[tf,V,^] e-^'fl. (2.12) 

JTxL s 

Here, the action is the sum S[U,if),iJj\ = S g [U] + S{[U,i/j,i/j] of the 0(a) improved 
plaquette action 

S*[U] = \'E l w(p)to(l-U(p)) (2.13) 
% p 

and the fermionic action 

&[tf,^=^^(z)(D + moMa:) (2.14) 

X 

for two degenerate flavours implicit in ip. For the special boundary conditions con- 
sidered below, the weight factor w(p) is the boundary improvement term c t [3] for 
time-like plaquettes at the boundary and one in all other cases. The value of c t has 
become available to 2-loop order in perturbation theory [2H] in the course of this 
work. Thus some of our data sets use the 1-loop value of c t , hence our simulations 
adopt two different actions. Because of universality we expect them to yield the same 
continuum limits for our observables. 

The 0(a) improved Wilson Dirac operator D includes the Sheikholeslami-Wohlert 
term [31] multiplied with the improvement coefficient c sw that has been determined 
with non-perturbative precision in |S2], and a boundary improvement term that is 
multiplied by the coefficient c t , which is known to 1-loop order. For details and 
notation we refer to |33|34j . 

The boundary conditions in the space directions are periodic for the gauge fields 
and periodic up to a global phase 9 for the fermion fields. The value of 9 was optimized 
at 1-loop order of perturbation theory [3U|. It turns out that a value close to 9 = n/5 
leads to a significantly smaller condition number of the fermion matrix than other 
values of 9 and thus to a smaller computational cost. Benchmarks in the relevant 
parameter range for our project have shown this too, and therefore we adopt the 
choice 9 = n/5 in this work. 

In the time direction, Dirichlet boundary conditions are imposed at Xq = and 
Xo = T. The quark fields at the boundaries are given by the Grassmann valued 
fields p, p at xo = and p', p' at xo = T, respectively. They are used as sources that 
are set to zero after differentiation. The gauge fields at the boundaries are chosen 
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such that a constant colour-electric background field, which is the unique (up to gauge 
transformations) configuration of least action, is generated in our space-time [I] . This 
is achieved by the diagonal colour matrices specified in ref . ^B] > parametrized by two 
dimensionless real parameters i] and v. 

A renormalized coupling g 2 may then be defined by differentiating the effective 
action T at the boundary point "A" of ^B] that corresponds to the choice rj = v = 0, 



dr] 



=S- ( 2 - 15 ) 



v =u=o 9 

The normalization k is chosen such that the tree-level value of g 2 equals g\ for all 
values of the lattice spacing. The boundary point "A" and especially the value v = 
are used, since the statistical error of the coupling turns out to be small for this 
choice. For general values of v we find another renormalized quantity v, 



dr 

dr] 



k { 1 - uv \ , (2.16) 



rj=0 
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that we have investigated as well to study the effects of dynamical fermions. 

The renormalized coupling depends on the system size, the lattice spacing a and 
on the quark mass. The bare quark mass mo is additively renormalized on the lattice 
because chiral symmetry is broken for Wilson fermions [SSI- Thus we define the bare 
mass by the PCAC relation that relates the axial current A^(x) = V ; ( a; )"y7At75' ? /'( a; ) to 
the pseudoscalar density P a (x) = ip(x) I ^'y5tp(x). Using the matrix elements /a and 
fp of A ^ and P a , respectively (cf. [33 ), the 0(a) improved PCAC mass is defined as 

/ x _ + $))/aQeq) + c A ad^d f P (x ) 
m[X0) ~ 2/p(z ) ■ [2A7) 

We have used 1-loop perturbation theory for c\ |34j. 

To fix all the details of our scheme, we define the bare current mass through 

( m(T/2), if T/a is even, 

\ [m((T - a)/2) + m((T + o)/2)] /2 , if T/a is odd. { ' 

This mass is tuned to zero, 

m(g ,m ,L/a) = 0, (2.19) 

so that we have a massless renormalization scheme, in which the only remaining 
external scale in the continuum limit is the system size L. 
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2.3 Computational Strategy 

In the last section we have defined the Schrodinger functional coupling g 2 . This 
finite- volume coupling runs with \x = 1 /L and - assuming monotonicity - there is a 
one-to-one relation between the value of the coupling and the system size L or energy 
scale jJL. In an abuse of notation we will from now on write g 2 (L) instead of g 2 (l/L). 

Our goal is to calculate the scale evolution of the strong coupling and the 
A-parameter of QCD in terms of a low-energy scale. We start the computation 
by choosing a value u for the renormalized coupling (which implicitly determines L) 
and by choosing a lattice resolution L/a. The theory can then be renormalized by 
tuning the bare parameters (3 = 6/g 2 and k — 1/(8 + 2am ) such that 

f(L)=u and m = 0. (2.20) 

Now we simulate a lattice with twice the linear size at the same bare parameters, 
that means at the same value of the lattice spacing, and thus with the physical extent 
2L, corresponding to a new renormalized coupling u' = g 2 (2L). 4 This determines the 
scale evolution of the renormalized coupling. It can be expressed through the lattice 
step scaling function 

Z(u,a/L)=?(2L)\ 32(L)=um=0 , (2.21) 
which is the key observable we compute. Finally, we obtain the step scaling function 

a(u) = lim £(u,a/L) (2.22) 

in the continuum limit by repeating these three steps with finer and finer lattice 
resolutions. This algorithm is iterated for a sequence of values for u to get the 
functional form of a(u). 

For small values of u the step scaling function cr(u) can be expanded in renor- 
malized perturbation theory, 

a(u) =u + s u 2 + Sim 3 + . • • , (2.23) 

4 The mass m on this larger lattice is different from zero by a lattice artifact which is expected 
to vanish in the continuum limit proportionally to a 2 . This has been verified in ref. |10j . 
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with the coefficients 



s = 2 6 ln2, (2.24) 
si = (2 6 hi2) 2 + 2 61 In 2, (2.25) 
s 2 = (2 6 ln2) 3 + 10 6 61 (In 2) 2 + 2 6 2 In 2. (2.26) 

The step scaling function cr(u) can be interpreted as an integrated discrete /3-function. 
Indeed, by using equation ()2.3|) we get 




(2.27) 



for the /3-function, which allows to calculate it recursively once the step scaling func- 
tion cr(-u) is known. 

To arrive at our main result, that is the A-parameter in terms of a low-energy 
scale, we solve the equation 

a (f(L/2)) = f{L) (2.28) 

recursively for g 2 (L/2). We start this recursion at a maximal value u max = g 2 (L maiX ) 
of the coupling. The value of u max is chosen such that the associated scale L max is 
a scale in the hadronic regime of QCD. Following the recursion (|2.28|) to larger and 
larger energies, we obtain the values for 

Ui = g 2 (2' l L max ) , i = 0,...,n, u = u max . (2.29) 

We perform n = 7 or n = 8 steps of this recursion and can in this way cover a 
scale separation of a factor 100 to 250. Eventually, for sufficiently large energies, 
perturbation theory can safely be applied. Then we use ()2.9|) with \i = 2 n /L max and 
with the /3-function truncated at 3-loop order, (|2.6|1 — (|2.8|) . The final result for AL max 
in the Schrodinger functional scheme can be converted to the MS scheme with (12.11)1 . 
We also check the admissibility of employing perturbation theory by studying the 
variation of our final result with respect to the number of non-perturbative steps n 
in the scale evolution of the strong coupling. 

2.4 Discretization effects 

The influence of the underlying space-time lattice on the evolution of the coupling can 
be estimated perturbatively j2Hj, by generalizing Symanzik's discussion |36|37|38j to 
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the present case. Close to the continuum limit we expect that the relative deviation 

$( u a /L) = ^ a l L )-< u ) = Sl ( a /L)u + 5 2 {a/L)u 2 + ... (2.30) 
a{u) 

of the lattice step scaling function from its continuum limit converges to zero with a 
rate roughly proportional to a/L. More precisely, since the action is 0(a) improved, 
we expect 

Si(a/L) ~ (d Q ,i + di,iln|) (|) 2 + ... , (2.31) 
5 2 {a/L) ~ e , 2 |+ ^o,2 + rfi, 2 ln| + rf 2 , 2 (ln|)^ (|) 2 + ... (2.32) 

for the 1-loop value of ct and the same form with eo,2 = for the 2-loop value of 
Ct. Note that the tree- level discretization effects vanish exactly, since we normalize 
the coupling such that its perturbative expansion starts with g$ for all values of the 
lattice spacing. 

The coefficients 5± and 5 2 are collected in table Q for the resolutions needed in 
this work. An expanded version of this table can be found in The entries in 



L/a 


Si 


rct=l-loop 
°2 


rct=2-loop 
°2 


4 


-0.0103 


0.0063 


-0.00007 


5 


-0.0065 


0.0049 


-0.00019 


6 


-0.0042 


0.0038 


-0.00041 


8 


-0.0021 


0.0029 


-0.00030 



Table 1: Discretization error of the step scaling function. 

the last column are very small. For larger values of L/a than shown in the table, 
^ct-2-ioop ^g^-ga^gg as expected. Since 52 t_1 ~ loop is of the order a/L, it is no surprise 
that it is much larger than 5 2 ~ 2 " loop - In fact, it is of the same size as 5\, for which 
the linear term in a/L is absent. 

The largest coupling at which the step scaling function has been computed with 
the 1-loop value of ct is u = 1.7319. With the 2-loop value of ct, this is u = 3.334. 
Tabled suggests that the step scaling function is only mildly affected by discretization 
effects. This will be demonstrated by our numerical results in section 0] 

We cancel the known perturbative cutoff effects for the respective actions by 
using 

S(2) K a l L ) = i W ^ /M 2 ( 2 - 33 ) 

v ' ' 1 + 8 1 (a/L)u + 5 2 {a/L)u 2 K ' 
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in the analysis of our Monte Carlo data. The perturbative estimate of the relative 
cutoff effects behaves as (a/L) x u 3 close to the continuum limit. 

2.5 Matching to a hadronic scheme 

As described so far, our computational strategy yields AL max , with L max defined by 
the value of the coupling itself. Since the latter is not experimentally measurable, it 
remains to relate L max to a hadronic scale. Here, the natural choice is the pion decay 
constant F n , since chiral perturbation theory provides an analytic understanding of 
the pion dynamics |40|41| , which is expected to help to control the extrapolations to 
the physical quark mass |H] as well as to infinite volume |42j . 

A computation of L max F % requires the knowledge of f{g$) = aF n (at a quark mass 
where m^/F^ takes its experimental value) and l(g$) = L max /a, where g 2 (L max ) = 
M ma X - We remind the reader that g 2 is defined at vanishing quark mass. The difference 
between the improved bare coupling g j^Sj and g is proportional to the light quark 
mass and can safely be neglected for physical values of the light quark mass. We 
therefore replace go ~ go in the following. The value of w max is restricted to be in 
the range covered by the computation of the scale dependence of the coupling and, 
for lattice spacings accessible in large- volume simulations, l(gl) should be sufficiently 
large. With both functions, I and /, defined for the same discretization, one finally 
wants to evaluate 

^max^| continuum = \M(gl)f(gl) . (2.34) 

Unfortunately, the results available in the literature |43|44j for f{g$) with our ac- 
tion [32] suffer from an uncertainty in the renormalization of the axial current, which 
has not yet been performed non-perturbatively. Also the 0(a) improvement of the 
current is known only perturbatively and it is not obvious that quark masses have 
been reached where chiral perturbation theory is applicable. 

At present, we thus prefer to relate L max to the frequently used hadronic radius r , 
which, according to phenomenological considerations, has a value of around 0.5 fm [13] 
and which has also been the reference scale in the zero-flavour theory, i.e. the quenched 
approximation [7j. Note that in a calculation in this approximation agreement was 
found (within its 3% precision) between F^r and the product of the experimental 
number for Fk times 0.5fm [Hj. Below, we thus evaluate L max /r = l(go)/p(go) where 
p(<7o) = r o/ a and translate to physical units via r = 0.5 fm. 
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3 Details of the numerical simulation and analysis 

In addition to the piece present in the pure gauge theory [TE], the central observ- 
able, dT/drj, receives a contribution due to the quark determinant. We describe its 
numerical evaluation by a stochastic estimator in appendix El Also detailed tables 
with simulation results are deferred to an appendix (JBj) . 



3.1 Tuning 

With a number of tuning runs we determine the bare parameters f3 and k such 
that eq. (|2.20|) is valid. Fulfilling the condition m(L) = precisely would require a 
fine tuning of the hopping parameter k. However, the uncertainty in £ owing to a 
small mismatch of m can be estimated perturbatively. To this end we compute the 
derivative of S with respect to z = mL, 



= <5>(a/L)u 2 + ... . (3.1) 

g 2 (L)=u,m(L)=z/L 



It turns out that $ is a slowly varying function of a/L. Thus, for our purpose it 
suffices to approximate it by its universal part, 



N t d 



0.00957iVf, (3.2) 



2=0 



where c\ t i(z) has been taken from [3U|. The typical precision in Monte Carlo sim- 
ulations is A(g~ 2 ) = g~ 4 A(g 2 ) = 0.003, as can be seen from the tables ITU1 and fTTI 
in appendix [BJ This means that an additional error of 0.001 x u 2 due to a slight 
mismatch of the mass m(L) is tolerable. Then it suffices to require 

am < 0.1-—. (3.3) 

The tables El and ^2 i n appendix El show that we have reached this precision in our 
simulations. 

Analogously, we stop the fine tuning of (3 if g 2 {L) = u well within the errors. 
Then we can correct for a small mismatch owing to g 2 {L) = u ^ u by using 

E(u, a/L) = E(u, a/L) + E'{u, a/L) x{u-u), (3.4) 



11 



with the perturbative estimate 

E'(u, a/L) = 9 ^ U : a/L) - ^ - 1 + 2.0- + Ss^ + W (3.5) 
ou ou 

for the derivative of the step scaling function E. This correction is always smaller 
than the statistical error of the step scaling function. 

Similarly, we convert the statistical error on u into an additional error of E, 

A (E(u)) « E'(u, a/L) x A(u) « (1 + 2s w + . . .) x A{u) . (3.6) 

This additional error is always much smaller than the error of g 2 (2L), to which it is 
added in quadrature. 

In some cases we have stopped the fine tuning of f3 and k after a number of runs 
and interpolated the results such that exactly the target coupling and mass zero with 
the errors as shown in the tables El and ^2 were obtained. 

The step scaling function E(w, a/L) and its partner with the perturbative cutoff 
effects being divided out, ^ 2 \u,a/L) (cf. ljT33l) ). are listed in table 121 



u 


L/a 


E(m, a/L) 


E( 2 )(u,a/L) 


u 


L/a 


E(u, a/L) 


E( 2 )(u,a/L) 


0.9793 


4 


1.0643(35) 


1.0686(35) 


1.5031 


4 


1.7204(56) 


1.7477(57) 




5 


1.0720(40) 


1.0738(41) 




5 


1.737(11) 


1.755(11) 




6 


1.0802(46) 


1.0807(46) 




6 


1.730(13) 


1.743(13) 




8 


1.0736(59) 


1.0729(59) 




8 


1.723(16) 


1.730(16) 


1.1814 


4 


1.3154(55) 


1.3199(56) 


2.0142 


4 


2.481(18) 


2.535(18) 




5 


1.3296(61) 


1.3307(61) 




5 


2.438(20) 


2.473(20) 




6 


1.3253(70) 


1.3249(70) 




6 


2.507(27) 


2.533(28) 




8 


1.3342(71) 


1.3323(71) 




8 


2.475(35) 


2.489(35) 


1.5031 


4 


1.7310(61) 


1.7332(61) 


2.4792 


4 


3.251(28) 


3.338(29) 




5 


1.756(12) 


1.754(12) 




5 


3.336(52) 


3.394(53) 




6 


1.745(12) 


1.741(12) 




6 


3.156(57) 


3.198(58) 












8 


3.326(52) 


3.351(53) 


1.7319 


4 


2.0583(76) 


2.0562(76) 


3.334 


4 


5.588(54) 


5.791(56) 




5 


2.083(21) 


2.076(21) 




5 


5.43(11) 


5.56(11) 




6 


2.058(20) 


2.049(20) 




6 


5.641(99) 


5.75(10) 












8 


5.48(13) 


5.53(13) 



Table 2: Step scaling functions S and E^ 2 '. The left hand side of the table contains the data with 
the 1-loop value of ct, while the data with c t = 2-loop are shown on the right. 
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3.2 Parameters 

Our choice of parameters is displayed in table El and table ^2 i n appendix [BJ The 
parameters shown are the results of a careful application of the tuning procedure 
explained in the last section. The tables reveal that the condition g 2 (L) = u is fulfilled 
to a good precision. The remaining deviations and the errors are then propagated 
into an additional error for as described above. 

3.3 Simulation costs and proper sampling of the configuration space 

Most of our results have been produced with the hybrid Monte Carlo algorithm 
(HMC) the polynomial hybrid Monte Carlo (PHMC) in the version proposed 
in |47|48| and for some of the L/a = 16 runs the hybrid Monte Carlo algorithm 
generalized to two pseudofermion fields (HMC 2 pf) |49|50| . 
We measure the cost of our simulations with the quantity 

-^cost = (update time in seconds on machine M) (3.7) 
x (error of l/g 2 ) 2 x (4a/T)(4a/L) 3 . 

In [51] the cost of a subset of the simulations discussed here (essentially up to L/a = 
12 and g 2 ~ 2.5) has been analyzed. A typical value for HMC at L/a = 10 and 
g 2 w 2.5 is M cost ~ 3.5, with respect to one board of APEmille [52] • To give an idea 
about the increase of the cost for the L/a = 16 simulations, we collect results from 
two different runs at the coupling g 2 « 2.5 in tableEJ We compute the autocorrelation 
time r int in units of trajectories in the way suggested in j5H]. 



algorithm 


step size 


Tint 


^update / [ S ] 


f 


iVtraj 


Accost 


HMC 


0.0625 


3.6(3) 


380 


2.46(5) 


4800 


6.9(6) 


HMC 2 pf 


0.111 


4.6(4) 


262 


2.55(5) 


5900 


5.6(5) 



Table 3: Cost estimates for two L/a = 16 runs; t U pdate/[s] is the time in seconds needed for one 
trajectory of length 1. The reference machine is an APEmille board. 

Related to the issue of estimating the autocorrelation time is always the question 
whether the algorithm samples the entire relevant configuration space efficiently If 
this is not the case, autocorrelation times may be largely underestimated and even 
systematically wrong results may be obtained. We now discuss two at least rough 
checks that our simulations do not suffer from such problems. 

1. To investigate the contributions to the coupling from sectors of configuration 
space with non-trivial topology, we have performed a set of simulations at the cou- 
pling g 2 « 2.5 with L/a = 8, 12. The topological charge Q(U) is determined through 
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a cooling procedure. Since sampling different topological sectors might be algorith- 
mically very difficult, we have employed the PHMC algorithm, whose flexibility can 
be exploited to enhance the transition rate among different sectors (at the price of 
increasing the fluctuations of the reweighting factor), and checked the results to be 
independent of the polynomial approximation used. 

Starting from a hot, random configuration, several of the independent replica 
were in a non-trivial topological sector. For these replica the smallest eigenvalue 
of D 2 (even-odd preconditioned) turns out to be one order of magnitude smaller 
than typical values in the topologically trivial sectors. After 0(100) trajectories for 
L/a = 8, respectively 0(1000) trajectories for L/a = 12, all the replica have zero 
topological charge and transitions to sectors with Q(U) ^ have not been observed 
in additional O(10 4 ) trajectories. 

From this we conclude that the PHMC algorithm can tunnel between different 
topological sectors, but for large L/a the transition rate is very small. In addition, 
the weight of the non-trivial sectors is too small for them to occur in a practical 
simulation at all (all tunnellings went to Q(U) = and none in the reverse direction). 
Their weight in the path integral is negligible. These statements have been checked 
for L/a = 8,12, and it appears safe to assume their validity also for larger L/a. 
Therefore, we decided to always start from a cold configuration, especially for the 
L/a = 16 simulations, to avoid thermalization problems. 

2. For the two largest couplings discussed here the distribution of dS/drj shows 
long tails toward negative values. The same effect was also observed in the computa- 
tion of the Schrodinger functional coupling in pure SU(3) gauge theory [T^j. We have 
related this tail to secondary local minima of the action [53] by measuring on cooled 
configurations the pure gauge contribution to the action Sg° o1 and to the coupling 
dSg° ol /dr]. This leads to metastabilities as shown for an L/a = 16 simulation with 
g 2 ~ 3.3 in figure ^ The upper panel is the Monte Carlo history (£md is the Monte 
Carlo time in units of molecular dynamics trajectories) of the gauge part of the ac- 
tion after cooling for two independent replica. The lower panel shows the history of 
dSg° ol /dr] for the same two replica. The correlation between metastable states and 
small (even negative) values of dSg° ol /dr] appears evident in this case. The action 
gcooi f Qr me tastability in the figure is consistent with the value for a secondary so- 
lution of the field equations given our choice for the boundary fields. Numerical 
evidence suggests that this solution is a local minimum. 

In order to estimate the weight of these contributions in our expectation values 
properly, we have enhanced their occurrence through a modified sampling similar 
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^ cool 




3000 



Figure 1: Monte Carlo histories of 5 1 ™ ' and dS^° ol /dr] from two independent replica (solid and 
dotted lines) in a simulation of a 16 4 lattice at g 2 ~ 3.3. 



to > adding to the HMC effective action a term 



7 



drj 



77=0 



+ — (7-7o) , 



(3.8) 



where 70 and are fixed to suitable (positive) values, while 7 is a dynamical variable. 
The expectation values in the original ensemble are then obtained by reweighting. 
By some deterministic cooling procedure with a fixed number of cooling steps we now 
define a quantity q whose value is 1 for metastable configurations, and otherwise, 
such that S q i + 8go = 1. For an arbitrary observable O we have an exact identity 



(O) = (5 ql O) + (S q0 O) 

= (5 ql O) + (0) I (l-(5 ql )) 



(3.9) 



with (O)i = (SgpO)/ (Sgo). If the main contribution to (O) comes from the config- 
urations with q = 0, a precise estimate of (O) just requires a precise estimate of 
(O)i, which can be obtained by an algorithm that samples only the q = sector, to- 
gether with rough estimates of (S q i) and {8 q \0) that can be obtained by the modified 
sampling. 

At g 2 ~ 3.3 we get (S q i) = 0.3(2)%, independent of L/a = 8, 12 within the error. 



15 



The effect of metastable states on the coupling is 0.10(2)%. At the coupling g 2 w 5.5 
the occurrence of metastable states is much more frequent, as expected, and we get 
(fiqi) — 5(1)%. At the same time their sampling is much easier already in the original 
ensemble, using either the PHMC or the HMC algorithm. In fact, we have repeated 
the L/a = 12 simulation for g 2 « 5.5 using PHMC for the ordinary ensemble (without 
f)3.8jl ) as in ^0], but measuring in addition the occurrence of metastable states. This 
turned out to be around 6%, and for g 2 we have obtained a result fully consistent 
with the number in [TU] . 

In summary, we can be confident that topologically non-trivial sectors are irrel- 
evant for our observables with our choice of parameters and at the present level of 
precision. In contrast, there are secondary minima in the action, which are visible as 
metastable states in the Monte Carlo sequence. They are relevant starting at g 2 ~ 3 
and have been taken into account efficiently by deviating from naive importance 
sampling and combining two different properly chosen ensembles. 

4 Results 

4.1 The strong coupling 

In figure 121 we show the approach of the step scaling function YP\u,a/L) to the 
continuum limit. The cutoff effects are small. Actually, all the data are compatible 
with constants. If we use simple fits to constants, the combined % 2 per degree of 
freedom for all the eight continuum extrapolations is about 1.4 regardless of the 
number of lattices included in the fit. Even the points at L/a = 4 are compatible 
with a constant continuum extrapolation. One possible strategy for the continuum 
extrapolation is thus a fit to a constant that uses the lattices with L/a = 6, 8. 

In this fit to constants we exclude the two coarsest lattices since there is always 
the danger of including systematic cutoff effects into the results coming from the 
lattices with large a/L. Therefore, we have also investigated two alternative fit pro- 
cedures. The first and most conservative one (denoted as "global fit" in the tables) 
is a combined continuum extrapolation of all the data sets, but excluding L/a = 4. 
Here we use the two-parameter ansatz 

S (2) (w,a/L) = a(u) +p x u 4 (a/L) 2 , X G {1-loop, 2-loop} , (4.1) 

for the lattice artifacts, where the coefficients p x are understood to be associated to 
the data with the 1-loop and the 2-loop value of c t , respectively. 5 This fit results 

5 We have also considered fits that use a common ansatz for the step scaling function for all the 
data sets. These lead to slightly smaller error bars of the final results. 
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in p l ~ loop = 0.08(13) and p 2 " loop = 0.01(4), quantifying that lattice artifacts are not 
detectable in our data. Moreover we have studied a mixed fit procedure, using a fit 
to constants for the lattices with L/a = 6,8 for the 2-loop improved data sets and 
the global fit ansatz 1)4.1)1 . which includes a slope for the cutoff effects, for the 1-loop 
improved data sets. 

These different fits are performed to investigate the uncertainties in the contin- 
uum results. All our plots refer to the most conservative of these three fit procedures, 
which leaves both p 1_loop and p 2_loop unconstrained. 



u 




a(u) 






global fit 


fit to constants 


mixed fit 


0.9793 


1.0736(44) 


1.0778(36) 


1.0736(44) 


1.1814 


1.3246(81) 


1.3285(50) 


1.3246(81) 


1.5031 


1.733(23) 


1.741(12) 


1.733(23) 


1.7319 


2.037(41) 


2.049(20) 


2.037(41) 


1.5031 


1.7440(97) 


1.738(10) 


1.738(10) 


2.0142 


2.488(26) 


2.516(22) 


2.516(22) 


2.4792 


3.311(52) 


3.281(39) 


3.281(39) 


3.3340 


5.60(16) 


5.670(80) 


5.670(80) 



Table 4: Continuum limit of the step scaling function. 



The results of the continuum limit extrapolation of the step scaling function 
YP\u,a/L) are recorded in table |U At u — 1.5031 we have two sets of data, one of 
which was produced with the 1-loop value and the other with the 2-loop value of ct. 
Both continuum results agree well within their errors, which is an independent check 
of our extrapolation procedures. 

We interpolate the values of tabled by a polynomial of degree 6 in u, the first 
coefficients up to u 3 being fixed by 2-loop perturbation theory, cf. ([2.23)1 . This inter- 
polation is depicted in figure El For small values of u < 2, the step scaling function is 
well described by perturbation theory. Actually, the perturbative lines shown in this 
plot are the solution of 

for cr(u) using the 2-loop and the 3-loop /3-function. Looking at the comparison of 
successive perturbative approximations to the non-perturbative results, it appears 
likely that higher orders would not improve the agreement at the largest coupling. 
Rather we appear to have reached a value for the coupling, where the perturbative 
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Figure 2: Continuum extrapolation of the step scaling function. 
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expansion has broken down. In fact, already in section 13.31 we have discussed the 
indications that fluctuations around a secondary minimum of the action are important 
at this value of the coupling. Such a mechanism may represent one possible source of 
non-perturbative effects. 6 

We adopt the parametrized form of the step scaling function to compute the 
A-parameter. To this end we start a recursion with a maximal coupling u max = 
(? 2 (L max ). For w max = 5.5, the recursive step (|2.28|) is solved numerically to get the 
couplings Ui, corresponding to the energy scales /i = 2 l /L max , that are quoted in 
table El We then insert these couplings into eq. ()2.9j) for the A-parameter, using 
there the 3-loop /^-function. This gives the results in the third column of table |S] 
Employing the 2-loop /3-function leads to results that are larger by roughly 0.02. The 
table shows that for u < 2 the A-parameter barely moves within its error bars. To 
be conservative, we use the global fit result and quote 

-ln(AL max ) = 1.09(7) at w max = 5.5 (4.3) 

as our final result, if the hadronic scale L max is defined through u max = 5.5. 

6 In this context we note further that for the pure gauge theory it has been demonstrated that 
the Schrodinger functional coupling grows exponentially with L at even larger values of L . We 
do not expect that any semiclassical picture is applicable in that regime but rather see this as a 
disorder phenomenon. 
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global fit constant fit, L/a = 6,8 mixed cont. extrap. 



X 






Ui 


— WAT" } 


Ui 


— lnf A T 





5.5 


0.957 


5.5 


0.957 


5.5 


0.957 


1 


3.306(40) 


1.071(25) 


3.291(18) 


1.081(12) 


3.291(19) 


1.081(12) 


2 


2.482(31) 


1.093(37) 


2.479(20) 


1.096(23) 


2.471(20) 


1.106(24) 


3 


2.010(27) 


1.093(48) 


2.009(19) 


1.096(35) 


2.003(19) 


1.106(35) 


4 


1.695(22) 


1.089(57) 


1.691(16) 


1.099(43) 


1.690(17) 


1.103(44) 


5 


1.468(18) 


1.087(65) 


1.462(14) 


1.109(49) 


1.464(15) 


1.100(52) 


6 


1.296(16) 


1.086(73) 


1.288(12) 


1.122(55) 


1.292(14) 


1.100(63) 


7 


1.160(14) 


1.086(82) 


1.151(11) 


1.138(62) 


1.157(13) 


1.101(74) 


8 


1.050(13) 


1.088(93) 


1.041(10) 


1.155(70) 


1.048(13) 


1.103(87) 



Table 5: Recursive computation of the A-parameter starting at uq = Umax = 5.5. 



We have in addition computed the A-parameter as a function of w max in the 
interval w max = 3.0 . . . 5.5. The results can be parametrized as 

- ln(AL max ) = — + ^ ln(6o w max ) - 0.1612 + 0.0379 « max . (4.4) 

M max ■zDq 

This parametrization is motivated by ()2.9j) and represents the central values of our 
data with a precision better than one permille. The absolute error of — ln(AL max ) 
that we quote for all the values of w max is 0.07. This means that we have calculated 
the A-parameter in units of L max with a precision of seven percent. 

The running of the Schrodinger functional coupling a(/x) = g 2 (1 / ft) / (47r) as a 
function of ///A is displayed in figure HJ The points refer to the entries of the second 
column of table El The symbol size is larger than their error. The difference between 
the perturbative and the non-perturbative running of the coupling looks small in this 
plot. However, if we had used perturbation theory only to evolve the coupling over 
the range considered here, the A-parameter would have been overestimated by up to 
14%, depending on w max . This corresponds to an extra error of 3% for the coupling in 
the range where its value is close to 0.12, corresponding to the physical value of o^jg 
at Mi- Needless to add, this error could of course not even be quantified without 
non-perturbative information. 

Furthermore, our non-perturbative coupling was designed to have a good per- 
turbative expansion, since we rely on the 3-loop /^-function in the (high-energy part 
of the) computation of the A-parameter. With our computation we have shown that 
for the Schrodinger functional coupling there is an overlapping region where both, 
perturbative and non-perturbative methods apply. In no way is this to be interpreted 
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Figure 4: Running of the strong coupling in the Schrodinger functional scheme. 
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as a general statement about QCD observables or couplings at certain energies. 

In figure|H]we show the non-perturbative /3-function in the Schrodinger functional 
scheme, together with the 2-loop and the 3-loop perturbation theory. It has been 
obtained recursively from ()2.27|) . The derivative of the step scaling function needed 
there has been calculated from the polynomial interpolating <j(u) (see continuous 
line in figure OJ). The non-perturbative data are fitted with two parameters beyond 
the 2-loop /^-function. The plot again shows an overlapping region in g, where the 
perturbative and the non-perturbative /5-functions agree well with each other. For 
a > 0.2, however, perturbation theory is no longer valid. Furthermore, the plot 
shows the difference between Nf = and Nf — 2. Already the leading coefficient bo 
of the /^-function depends on the number of flavours, and this is nicely reflected in 
the figure. 

4.2 Computation of v as a function of the strong coupling 

The difference between the quenched approximation and the two-flavour theory is 
also apparent in the renormalized quantity v defined in f)2 . 1 6j) . As a function of the 
coupling u we write (at zero quark mass) 

v = u(u) = lim Q(u,a/L). (4-5) 

a/L—>0 

In perturbation theory, Q is known to 2-loop order, 

tt(u, a/L) = (vi + v 2 u) (1 + e 1 (a/L) + e 2 (a/L)u) + 0(u 2 ) . (4.6) 
Here |io'K()l^| . 

vi = 0.0694603(1) + 0.0245370(1)^, (4.7) 
v 2 = -0.001364(14) - 0.000101(17)iV f - 0.0003362(30)iV f 2 , (4.8) 

and the perturbative cutoff effects are listed in table |H1 Note that the tree- level coeffi- 
cient of v vanishes exactly because of the definition of the couplings. The perturbative 
results indicate a large effect for going from the zero- to the two-flavour theory. 

The first step in the analysis is to project v on zero mass. To this end we 
have obtained the crude estimate dv/d(am) m —0.15(4) at constant u from the 
matching runs at the smallest coupling and with L/a = 4. We use it also at the 
other couplings. After this projection, the perturbative cutoff effects are eliminated 
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= 


N f = 


= 2 


L/a 


ei(a/L) 


e 2 (a/L) 


ei(a/L) 


e 2 (a/L) 


4 


0.1880609(17) 


-0.02020(24) 


0.3044344(25) 


0.03725(43) 


5 


0.1085045(16) 


-0.01347(22) 


0.1822083(22) 


0.01255(39) 


6 


0.0677292(15) 


-0.00910(22) 


0.1114033(21) 


0.00325(36) 


7 


u.u4ouyyo^io j 


n c\c\p.p.c\( r >^ ^ 

— U.UUDDU^Zl ) 


u.u i Z i oDO^ZU ) 


— U.UUU4 i ) 


8 


0.0335967(15) 


-0.00513(21) 


0.0502340(20) 


-0.00130(34) 


10 


0.0203927(15) 


-0.00357(21) 


0.0280832(19) 


-0.00118(34) 


12 


0.0138002(15) 


-0.00273(20) 


0.0181305(19) 


-0.00091(33) 


14 


0.0099904(15) 


-0.00219(20) 


0.0127946(19) 


-0.00074(33) 


16 


0.0075781(15) 


-0.00181(20) 


0.0095635(19) 


-0.00064(33) 


20 


0.0047983(14) 


-0.00133(20) 


0.0059634(19) 


-0.00050(33) 


24 


0.0033132(14) 


-0.00103(20) 


0.0040873(19) 


-0.00040(33) 




Table 6: 


Cutoff effects of v in 


perturbation theory. 





(similarly to ()2.33|) ) by replacing Q by 



Q(u, a/L) 



1 + ei (a/L) +e 2 (a/L)u 



(4.9) 



in the analysis. In contrast to the coupling, this correction is substantial and the 
resulting continuum extrapolation is much smoother |10j . 

Then we project QS 2 \u,a/ L) on some reference couplings in the range u = 
0.9793 . . . 5.5, using a numerical estimate for the slope. In principle, we would have 
to propagate the error of u = g 2 into an extra error of il^(u, a/L). However, it turns 
out, both from perturbation theory and from the non-perturbative fits later, that the 
variation of v with u is so small that this extra error can safely be neglected. 

We make an ansatz linear in (a/L) 2 for the continuum extrapolation. The data 
for the two different actions are extrapolated separately and the continuum results 
are then averaged according to their weight. The data at L/a = 4, 5 are left out. 

The results are displayed in figure El together with v at N{ = from ref. |16j . 
The comparison shows that v increases by almost a factor two when going from the 
quenched approximation to Nf = 2. 

4.3 Evaluation of Ar 

For the O(a) improved action (321 used in our computations, the low-energy scale r 
has been calculated at (3 = 6/^q = 5.2 (or a « 0.1 fm) by two groups |56|44j . Recently, 
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Figure 6: v as a function of the strong coupling a. The dashed lines show the 2-loop perturbation 
theory while the solid lines are non-perturbative fits for a < 0.3 including an additional a 2 -term. 

these large-volume Nf = 2 simulations of QCD have been extended to smaller values 
of the lattice spacing, namely 5.2 < (3 < 5.4 |2Zj. In order to obtain Aro and to 
study its dependence on the lattice spacing, we here use results for r Q /a of [HZ] at 
(3 = 5.2, 5.29, 5.4 and compare also to the numbers resulting from r /a of [44J. 

First, we obtain the renormalized coupling on lattices with extent L/a = 4, 6, 8 
at the three chosen values of (3. It is listed in table [7| The hopping parameters k 
are taken from |5T. They correspond to roughly massless pions and thus massless 
quarks. We checked that reasonable changes of k, e.g. requiring eq. (|2.19|) . affect our 
analysis only to a negligible amount. We then set the improvement coefficient c t to 
its 2-loop value and obtain L mBX / a for the three values of (3 combined with two fixed 
values = ^ 2 (-^max) = 3.65 and M max = 4.61 by an interpolation of the data in 
table [7| These values of L max / a, which are recorded in table |Hl are insensitive to the 
interpolation formula used. 

Second, we analyze the raw data for r^/a at finite bare quark masses, m q , in order 
to obtain the value corresponding to massless quarks (the up- and down-quark masses 
may safely be neglected in this context). At each bare coupling, three different quark 
masses, corresponding to pion masses between about 500 MeV and 1 GeV, have been 
simulated in [57| . As seen in figure [7[ the radius r depends approximately linearly 
on the bare quark mass in this range. The figure also demonstrates that the slope 
is strongly cutoff dependent; its magnitude decreases quite rapidly as (3 increases 
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Ct=l-loop c t =2-loop 



p 




K 


L/a 


9\L) 


5.20 


0. 


13600 


4 


3.32(2) 


3.65(3) 


5.20 


0. 


13600 


6 


4.31(4) 


4.61(4) 


5.29 


0. 


13641 


4 


3.184(16) 


3.394(17) 


5.29 


0. 


13641 


6 


4.059(32) 


4.279(37) 


5.29 


0. 


13641 


8 


5.34(8) 


5.65(9) 


5.40 


0. 


.13669 


4 


3.016(20) 


3.188(24) 


5.40 


0. 


13669 


6 


3.708(31) 


3.861(34) 


5.40 


0. 


13669 


8 


4.704(59) 


4.747(63) 



Table 7: Simulation results for g 2 (L) at low (3. The hopping parameters k are set to the critical 
ones (k c ) of |57) . 







"max 


= 3.65 


"max 


= 4.61 


p 


r /a 


Lm&x/ ®> 




Lxaax./ ®> 




5.20 
5.29 
5.40 


5.45(5)(20) 
6.01(4)(22) 
7.01(5)(15) 


4.00(6) 
4.67(6) 
5.43(9) 


0.655(27) 
0.619(25) 
0.621(17) 


6.00(8) 
6.57(6) 
7.73(10) 


0.610(25) 
0.614(24) 
0.609(16) 



Table 8: Low-energy scale tq in the chiral limit and the combination A^jg r$ as obtained for two 
values of u max = .g 2 (L max ). 

(the lattice spacing decreases). 7 This has been noted earlier [HE] and reminds us 
that the study of lattice artifacts is important. On the other hand, we are here not 
interested in the slope but in r at zero quark mass. We estimate it by a simple 
linear extrapolation in m q . Since the linear behaviour is not guaranteed to extend 
all the way to zero quark mass, we include a systematic error of the extrapolation in 
addition to the statistical one: an uncertainty of half the difference of the last data 
point and the extrapolated value is added as the second error in table |H1 Within the 
total error, dominated by the systematic one due to the extrapolation, our values of 

7 The slope visible in figure is not directly a physical observable, since (up to a-effects) the 
renormalized quark mass is given by tor = Z m m q rather than by m q . However, it appears unlikely 
that the strong dependence of the slope on go is cancelled by Z m , since the latter is expected to 
be a weak function of the bare coupling. More details can be found in |58j . where also a properly 
renormalized slope has been analyzed. 
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Figure 7: Dependence of ro on the bare quark mass m q = (1/k— 1/k c )/2. Both quantities are 
rescaled (made dimensionless) by the extrapolated value of ro, denoted by ro(0). The uncertainty of 
this rescaling is not propagated into the errors. For k c , we take the values of [SZ], listed in tabled 



r /a do agree with those quoted in |56|44p57| . where somewhat different ansatze - 
and in 03] also different data - have been used. 

Third, combining eqs. ()4.4|) and (j2.11|) with r /a and L mSuX / a of table |HJ one 
arrives at the columns A^jg r in that table. Here we have added errors in quadrature, 
except for the one attributed to eq. ()4.4|) . which is independent of the bare coupling. 
The resulting Aj^g ro are remarkably stable with respect to the change of the lattice 
spacing and also the choice of M max - In particular, for L max /a > 4.5 all numbers 
including their errors are covered by the interval 0.58 < Aj^gr < 0.66 and also the 
central value obtained with the worst discretization L max /a = 4 is inside. We thus 
quote 

A^r = 0.62(4) (4) (4.10) 

as our result, the second error deriving from the 7% error on AL max . 

We finally mention that we repeated the above analysis also for c t to 1-loop 
precision and the values M max = 3.32 and 4.31 suggested by table The central 
values for A^ ro are then up to 15% lower than for c t at 2-loop precision, but this 
difference shrinks as L max /a grows. At the largest value of L max /a, the number for 
Af^g r is again fully contained in the range A^jg r = 0.62(4). 
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5 Conclusions 



This non-perturbative QCD computation has required extensive simulations of Nf = 2 
QCD with 0(a) improved massless Wilson fermions in finite volume. In our situation, 
discretization errors turned out to be very small, as seen for example in figure [2] and 
also in [!(]]. Although we could achieve the necessary precision only on lattices up to 
a size 16 4 , the smallness of discretization errors allowed to obtain the running of the 
QCD coupling in the continuum limit and to good accuracy. As in the Nf — (pure 
gauge) theory, the energy dependence of the coupling in the Schrodinger functional 
scheme is now known over more than two orders of magnitude in the energy scale. 
This is the main result of our investigation. 

The iVf-dependence is best illustrated in figure where we also observe excellent 
agreement with 3-loop perturbation theory for a < 0.2, while for larger a, the non- 
perturbative /3-function breaks away rapidly from the perturbative approximation. At 
a ~ 0.3 . . . 0.4, a couple of additional higher-order perturbative terms with coefficients 
of a reasonable size would not be able to come close to the non-perturbative j3- 
function. 

To calibrate the overall energy scale, one fixes a large enough value of the coupling 
to be in the low-energy region and relates the associated distance, L max , to a non- 
perturbative, large- volume observable. For technical reasons, explained in section 12*31 
we have chosen the hadronic radius r , which has an unambiguous definition in terms 
of the force F(r) between static quarks, via rgF(r ) = 1.65 [""""]. This quantity has 
also been chosen in the computation of the A-parameter for Nf = [7j. We compare 
Aro for various numbers of flavours in table [21 



N t 


A MS r ° 


reference 




remarks 







0.60(5) 


7j 
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0.62(4)(4) 


this work 
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0.57(8) 


|5f)l60l61j 


DIS 


NNLO & r = 


= 0.5 fm 


4 


0.74(10) 


□ 


world 


average & r 


= 0.5 fm 


5 


0.54(8) 


M 


world 


average & r 


= 0.5 fm 



Table 9: The QCD A-parameter in units of rg. Non-perturbative, purely theoretical determinations 
for Nf = 0, 2 are compared to extractions of A from high-energy scattering experiments, using high- 
order perturbation theory combined with the phenomenological estimate vq w 0.5 fm |45| . 

The last two entries in the table represent one and the same world average by 
S. Bethke of a- measurements. They are related by the perturbative matching of the 
effective theories with Nf = 4 and Nf = 5 massless quarks ■ While little can be said 
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on the iVf-dependence of A™ r on general grounds, the prediction that there should 
be a significant drop from Nf = 4 to AT f = 5 depends only on perturbation theory at 
the scale of the b-quark mass and should thus be reliable. A similar statement for 
the change from Nf = 3 to Nf = 4 is less certain as it involves physics at and below 
the mass of the charm quark. 

Another relevant issue in the above comparison is that we used tq « 0.5 fm to 
relate the high-energy experiments to our theoretical predictions. Although it appears 
unlikely that r differs by 10% from this value, a true error is difficult to estimate until 
a reliable non-perturbative computation of e.g. r F n has been performed. Indeed, 
such a computation, or more directly the computation of L max x F n for Nf — 2, is 
the most urgent next step to be taken in our programme. After that, the effect of 
the remaining (massive) quarks needs to be estimated. 

Keeping the above caveats in mind, we still may convert the A-parameter to 
physical units and obtain 

= 245(16)(16)MeV [with r = 0.5 fin] ■ ( 5A ) 

Although in this case the four-flavour theory has not yet been reached, it is a very 
non-trivial test of QCD that the non-perturbative results, which use experimental 
input at low energies of order 1/Yo ~ 400 MeV, agree roughly with the high-energy, 
perturbative extractions of A. Unravelling the details in this comparison will still 
require some work; some of it was just mentioned. 

Now, that cx(n) is known, the tables presented in this work also provide the bare 
parameters of our lattice action needed in the computation of the energy dependence 
of the renormalized quark mass and composite operators. These are then readily 
related to the appropriate renormalization group invariants. 

Acknowledgements 

We are indebted to Martin Luscher who founded the ALPHA Collaboration and who 
led ground-breaking work for this project - as demonstrated by the references we 
quote. We further thank Achim Bode, Bernd Gehrmann, Martin Hasenbusch, Karl 
Jansen, Francesco Knechtli, Stefan Kurth, Hubert Simma, Stefan Sint, Peter Weisz 
and Hartmut Wittig for many useful discussions and collaboration in early parts of 
this project JUj- We are grateful to Gerrit Schierholz for communicating the results 
of |57] • We further thank DESY for computing resources and the APE Collaboration 
and the staff of the computer centre at DESY Zeuthen for their support. 

The computation of a s is one project of SFB Transregio 9 "Computational Par- 
ticle Physics" and has been strongly supported there as well as in Graduiertenkolleg 



28 



GK 271 by the Deutsche Forschungsgemeinschaft (DFG). We thank our colleagues 
in the SFB for discussions, in particular Johannes Blumlein, Kostia Chetyrkin and 
Fred Jegerlehner. This work has also been supported by the European Community's 
Human Potential Programme under contract HPRN-CT-2000-00145. 



A Evaluation of dT/drj 

Our central observable, eq. (|2.15|) . translates into the expectation value 



£ = (£) = /^) + (*f) , (A.i) 

where the pure gauge part dSJdr] has been discussed in [16^ and 

0-77 d?7 \ dr] / v 

Here we have used Sf s = — A/fTrlnQ with Q = 75 (-D + mo) = and D the (one- 
flavour) Dirac operator including improvement terms. As usual, Sf ff is obtained after 
integrating out the fermion fields and the trace extends over colour and Dirac indices 
as well as over the space-time points. The last expression in eq. flA.2|) represents an 
average over a complex random field <p(x) with the property 

( ( Ptx( x ) ( Pd$(v)) v = $cd <W S xy , (A.3) 

where c,d denote colour indices and a,/3 are Dirac indices. We may finally rewrite 
dSf /dr] in the form 

% = - N <( E itr!i«|^[^)x» W + x W ^W] 

x\xo£{a,T— a} 

V, (A.4) 

where we have used the fact that dQ/drj vanishes except for the clover terms, which 
are diagonal in the coordinate x and contribute only on the time slices Xq = a and 
Xq = T — a. Now "tr" is over spin and colour only. Eq. (jA.4|) is in the form used in 
[63J for the contribution of the clover terms to the pseudofermionic force in the HMC 
algorithm and is evaluated analogously. Only one solution of the Dirac equation is 
needed. 

Of course, the average over the gauge fields can be interchanged with the one 
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over the field (p and one may replace (|A.2|) by an average over a finite number of 
fields ip drawn from a distribution satisfying ljA.3jl . We found that the fluctuations 
of such a noisy (unbiased) estimator for dSf s /dr) are small compared to the ones of 
dS s /dr}, already when only one field (p (from a Gaussian distribution) is used per 
gauge configuration. This has hence been our method of choice in all simulations. 

B Detailed numerical results 

The following tables list detailed parameters and results of our simulations. 



30 



L/a (3 k g 2 A(g 2 ) v A(v) m A(m) 

u = 0.9793 

4 9.2364 0.1317486 0.9793 0.0007 0.1557 0.0019 -0.00600 O.OOOTI 

5 9.3884 0.1315391 0.9794 0.0009 0.1322 0.0023 0.00197 0.00005 

6 9.5000 0.1315322 0.9793 0.0011 0.1266 0.0016 -0.00014 0.00003 
8 9.7341 0.131305 0.9807 0.0017 0.1177 0.0042 0.00074 0.00006 
8 9.2364 0.1317486 1.0643 0.0034 0.1244 0.0061 0.00010 0.00004 
10 9.3884 0.1315391 1.0721 0.0039 0.1151 0.0077 0.00210 0.00003 
12 9.5000 0.1315322 1.0802 0.0044 0.1227 0.0072 -0.00091 0.00002 
16 9.7341 0.131305 1.0753 0.0055 0.1047 0.0080 -0.00008 0.00003 
u = 1.1814 

4 8.2373 0.1327957 1.1814 0.0005 0.1483 0.0016 0.00100 0.00011 

5 8.3900 0.1325800 1.1807 0.0012 0.1353 0.0018 -0.00018 0.00009 

6 8.5000 0.1325094 1.1814 0.0015 0.1269 0.0014 -0.00036 0.00003 
8 8.7223 0.1322907 1.1818 0.0029 0.1141 0.0048 -0.00115 0.00004 
8 8.2373 0.1327957 1.3154 0.0055 0.1209 0.0061 0.00020 0.00005 
10 8.3900 0.1325800 1.3287 0.0059 0.1128 0.0070 0.00097 0.00007 
12 8.5000 0.1325094 1.3253 0.0067 0.1304 0.0068 -0.00102 0.00002 
16 8.7223 0.1322907 1.3347 0.0061 0.1065 0.0049 -0.00194 0.00002 
u = 1.5031 

4 7.2103 0.1339411 1.5031 0.0010 0.1437 0.0010 -0.00074 0.00010 

5 7.3619 0.1339100 1.5044 0.0027 0.1250 0.0031 0.00052 0.00010 

6 7.5000 0.1338150 1.5031 0.0025 0.1201 0.0024 -0.00078 0.00004 
8 7.2103 0.1339411 1.7310 0.0059 0.1151 0.0037 0.00959 0.00004 
10 7.3619 0.1339100 1.7581 0.0113 0.1062 0.0084 0.00257 0.00005 
12 7.5000 0.1338150 1.7449 0.0119 0.1223 0.0073 -0.00138 0.00004 

u = 1.7319 

4 6.7251 0.1347424 1.7319 0.0020 0.1378 0.0009 -0.00181 0.00013 

5 6.8770 0.1346900 1.7333 0.0032 0.1272 0.0025 -0.00005 0.00011 

6 7.0000 0.1345794 1.7319 0.0034 0.1161 0.0023 -0.00002 0.00005 
8 677251 0.1347424 2.0583 0.0070 0.1008 0.0032 0.01051 0.00005 
10 6.8770 0.1346900 2.0855 0.0208 0.0934 0.0080 0.00335 0.00006 
12 7.0000 0.1345794 2.0575 0.0196 0.0833 0.0082 -0.00049 0.00006 



Table 10: Simulation parameters and results using the 1-loop value for c t . 
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L/a (3 k g 2 AQ7 2 ) v A(v) m A(m) 

u = 1.5031 

4 7.2811 0.1338383 1.5031 0.0012 0.1434 0.0013 0.00043 0.00015 

5 7.4137 0.1338750 1.5033 0.0026 0.1310 0.0027 -0.00083 0.00009 

6 7.5457 0.1337050 1.5031 0.0030 0.1247 0.0031 0.00072 0.00008 
8 7.7270 0.133488 1.5031 0.0035 0.1219 0.0032 -0.00084 0.00003 
8 7.2811 0.1338383 1.7204 0.0054 0.1091 0.0037 0.00928 0.00004 
10 7.4137 0.1338750 1.7372 0.0104 0.1110 0.0063 0.00109 0.00004 
12 7.5457 0.1337050 1.7305 0.0122 0.0893 0.0079 -0.00006 0.00008 
16 7.7270 0.133488 1.7231 0.0151 0.1017 0.0092 -0.00154 0.00019 
u = 2.0142 

4 6.3650 0.1353200 2.0142 0.0024 0.1349 0.0017 0.00000 0.00023 

5 6.5000 0.1353570 2.0142 0.0044 0.1236 0.0024 0.00002 0.00011 

6 6.6085 0.1352600 2.0146 0.0056 0.1205 0.0029 0.00030 0.00009 
8 6.8217 0.134891 2.0142 0.0102 0.0991 0.0045 0.00049 0.00007 
8 6.3650 0.1353200 2.4814 0.0172 0.1016 0.0049 0.01318 0.00008 
10 6.5000 0.1353570 2.4383 0.0188 0.0900 0.0050 0.00367 0.00005 
12 6.6085 0.1352600 2.5077 0.0259 0.1074 0.0069 0.00013 0.00004 
16 6.8217 0.134891 2.475 0.031 0.0916 0.0076 -0.00053 0.00006 
u = 2.4792 

4 5.8724 0.1360000 2.4792 0.0034 0.1206 0.0016 0.00000 0.00026 

5 6.0000 0.1361820 2.4792 0.0073 0.1085 0.0023 -0.00009 0.00014 

6 6.1355 0.1361050 2.4792 0.0082 0.1025 0.0032 0.00000 0.00013 
8 6.3229 0.1357673 2.4792 0.0128 0.1015 0.0053 0.00000 0.00016 
8 5.8724 0.1360000 3.2511 0.0277 0.0859 0.0043 0.01819 0.00011 
10 6.0000 0.1361820 3.3356 0.0502 0.0796 0.0064 0.00579 0.00009 
12 6.1355 0.1361050 3.1558 0.0552 0.0801 0.0079 0.00078 0.00007 
16 6.3229 0.1357673 3.3263 0.0472 0.0806 0.0074 0.00039 0.00017 
u = 3.3340 

4 5.3574 0.1356400 3.3340 0.0109 0.1087 0.0013 0.00000 0.00040 

5 5.5000 0.1364220 3.3340 0.0182 0.0965 0.0018 -0.00004 0.00017 

6 5.6215 0.1366650 3.3263 0.0196 0.0894 0.0031 0.00051 0.00018 
8 5.8097 0.1366077 3.334 0.019 0.0887 0.0042 0.00000 0.00004 
8 5.3574 0.1356400 5.588 0.049 0.0576 0.0036 0.03163 0.00019 
10 5.5000 0.1364220 5.430 0.098 0.0585 0.0048 0.01126 0.00012 
12 5.6215 0.1366650 5.624 0.089 0.0607 0.0043 0.00334 0.00015 
16 5.8097 0.1366077 5.4763 0.1236 0.0689 0.0052 0.00048 0.00004 



Table 11: Simulation parameters and results using the 2-loop value for c t . 



32 



References 

[1] S. Bethke, Nucl. Phys. Proc. Suppl. 135 (2004) 345, hep-ex/0407021. 

[2] K.I. Ishikawa, plenary talk at LATTICE '04, Fermilab, Batavia IL, USA, 2004, 
hep-lat/0410050. 

[3] M. Liischer, P. Weisz and U. Wolff, Nucl. Phys. B359 (1991) 221. 

[4] M. Liischer et al, Nucl. Phys. B384 (1992) 168, hep-lat/9207009. 

[5] S. Sint, Nucl. Phys. B451 (1995) 416, hep-lat/9504005. 

[6] ALPHA, G. de Divitiis et al., Nucl. Phys. B437 (1995) 447, hep-lat/9411017. 

[7] ALPHA, S. Capitani et al, Nucl. Phys. B544 (1999) 669, hep-lat/9810063. 

[8] ALPHA, J. Garden et al, Nucl. Phys. B571 (2000) 237, hep-lat/9906013. 

[9] ALPHA, J. Heitger et al., Nucl. Phys. Proc. Suppl. 106 (2002) 859, hep- 
lat/0110201. 

[10] ALPHA, A. Bode et al., Phys. Lett. B515 (2001) 49, hep-lat/0105003. 

[11] TXL, A. Spitz et al., Phys. Rev. D60 (1999) 074502, hep-lat/9906009. 

[12] QCDSF, M. G6ckeler et al., (2004), hep-lat/0409166. 

[13] C. Davies et al., Nucl. Phys. Proc. Suppl. 119 (2003) 595, hep-lat/0209122. 

[14] B. Bunk et al., Nucl. Phys. B697 (2004) 343, hep-lat/0403022. 

[15] M. Liischer et al., Nucl. Phys. B389 (1993) 247, hep-lat/9207010. 

[16] M. Liischer et al., Nucl. Phys. B413 (1994) 481, hep-lat/9309005. 

[17] S. Sint, Nucl. Phys. Proc. Suppl. 42 (1995) 835, hep-lat/9411063. 

[18] R. Sommer, Non-perturbative Renormalization of QCD, Lectures given at 36th 
Internationale Universitatswochen fur Kernphysik und Teilchenphysik, Schlad- 
ming, Austria, 1997, hep-ph/9711243. 

[19] M. Liischer, Advanced lattice QCD, Lectures given at Les Houches Summer 
School in Theoretical Physics, Les Houches, France, 1997, hep-lat/9802029. 

[20] C.G. Callan Jr., Phys. Rev. D2 (1970) 1541. 

33 



[21] K. Symanzik, Commun. Math. Phys. 18 (1970) 227. 

[22] S. Weinberg, Phys. Rev. D8 (1973) 3497. 

[23] Y. Frishman, R. Horsley and U. Wolff, Nucl. Phys. B183 (1981) 509. 

[24] O.V. Tarasov, A. A. Vladimirov and A.Y. Zharkov, Phys. Lett. B93 (1980) 429. 

[25] M. Liischer and P. Weisz, Nucl. Phys. B452 (1995) 213, hep-lat/9504006. 

[26] M. Liischer and P. Weisz, Nucl. Phys. B452 (1995) 234, hep-lat/9505011. 

[27] B. Alles, A. Feo and H. Panagopoulos, Nucl. Phys. B491 (1997) 498, hep- 
lat/9609025. 

[28] A. Bode and H. Panagopoulos, Nucl. Phys. B625 (2002) 198, hep-lat/0 110211. 

[29] ALPHA, A. Bode, P. Weisz and U. Wolff, Nucl. Phys. B576 (2000) 517, Erratum- 
ibid. B600 (2001) 453, Erratum-ibid. B608 (2001) 481, hep-lat/9911018. 

[30] S. Sint and R. Sommer, Nucl. Phys. B465 (1996) 71, hep-lat/9508012. 

[31] B. Sheikholeslami and R. Wohlert, Nucl. Phys. B259 (1985) 572. 

[32] ALPHA, K. Jansen and R. Sommer, Nucl. Phys. B530 (1998) 185, hep- 
lat/9803017. 

[33] M. Liischer et al, Nucl. Phys. B478 (1996) 365, hep-lat/9605038. 

[34] M. Liischer and P. Weisz, Nucl. Phys. B479 (1996) 429, hep-lat/9606016. 

[35] K.G. Wilson, Phys. Rev. D10 (1974) 2445. 

[36] K. Symanzik, Some topics in quantum field theory, in Mathematical problems 
in theoretical physics, eds. R. Schrader et al., Lecture Notes in Physics Vol. 153 
(Springer, New York, 1982). 

[37] K. Symanzik, Nucl. Phys. B226 (1983) 187. 

[38] K. Symanzik, Nucl. Phys. B226 (1983) 205. 

[39] B. Gehrmann et al., Nucl. Phys. B612 (2001) 3, hep-lat/0106025. 

[40] S. Weinberg, Physica A96 (1979) 327. 

[41] J. Gasser and H. Leutwyler, Ann. Phys. 158 (1984) 142. 



34 



[42] G. Colangelo, plenary talk at LATTICE '04, Fermilab, Batavia IL, USA, 2004, 
hep-lat/0409111. 

[43] UKQCD, A.C. Irving et al., Phys. Lett. B518 (2001) 243, hep-lat/0107023. 

[44] JLQCD, S. Aoki et al., Phys. Rev. D68 (2003) 054502, hep-lat/0212039. 

[45] R. Sommer, Nucl. Phys. B411 (1994) 839, hep-lat/9310022. 

[46] S. Duane et al, Phys. Lett. B195 (1987) 216. 

[47] R. Frezzotti and K. Jansen, Nucl. Phys. B555 (1999) 395, hep-lat/9808011. 

[48] R. Frezzotti and K. Jansen, Nucl. Phys. B555 (1999) 432, hep-lat/9808038. 

[49] M. Hasenbusch, Phys. Lett. B519 (2001) 177, hep-lat/0107019. 

[50] M. Hasenbusch and K. Jansen, Nucl. Phys. B659 (2003) 299, hep-lat/0211042. 

[51] ALPHA, R. Frezzotti et al, Comput. Phys. Commun. 136 (2001) 1, hep- 
lat/0009027. 

[52] A. Bartoloni et al., Nucl. Phys. Proc. Suppl. 106 (2002) 1043, hep-lat/01 10153. 

[53] ALPHA, U. Wolff, Comput. Phys. Commun. 156 (2004) 143, hep-lat/0306017. 

[54] ALPHA, M. Delia Morte et al., Nucl. Phys. Proc. Suppl. 119 (2003) 439, hep- 
lat/0209023. 

[55] A. Bode, v to two loop including cutoff effects, unpublished notes, 2001. 
[56] UKQCD, C.R. Allton et al., Phys. Rev. D65 (2002) 054502, hep-lat/0107021. 
[57] QCDSF, M. Gockeler et al., (2004), hep-ph/0409312. 

[58] ALPHA, R. Sommer et al., Nucl. Phys. Proc. Suppl. 129 (2004) 405, hep- 
lat/0309171. 

[59] J. Blumlein, H. BSttcher and A. Guffanti, (2004), hep-ph/0407089. 

[60] S. Moch, J.A.M. Vermaseren and A. Vogt, Nucl. Phys. B688 (2004) 101, hep- 
ph/0403192. 

[61] W.L. van Neerven and E.B. Zijlstra, Phys. Lett. B272 (1991) 127. 

[62] W. Bernreuther and W. Wetzel, Nucl. Phys. B197 (1982) 228. 

[63] K. Jansen and C. Liu, Comput. Phys. Commun. 99 (1997) 221, hep-lat/9603008. 



35 



